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IMPLICIT  SCHEME  FOR  PLASMA  TRANSPORT  IN  3D 
S.  Dinkevich*  and  A.  Bayliss** 
1.  Detailed  formulation  of  anisotropic  plasma  transport  in  3D  is 
given  in  [1,2,3].  All  notation  is  taken  from  there.   Plasma  transport 
in  3D  is  governed  by  the  following  three  parabolic  equations: 


9y 
3t 


av   3v2 
F(>jJ^  ,v,  -f-,  -^]    ,  (1) 


i|j^=const        ~   3^^   3^2' 


where 


V  =  (v^  ,Ci ,0q)^  .  (2) 

Normalization  is  defined  below: 

IjJ-  is  the  normalized  toroidal  flux, 

v-j  is  the  normalized  rotational  transform, 

Ci  is  the  normalized  relative  plasma  density, 

0  is  the  normalized  entropy  density, 

t  is  non-dimensional  time. 

This  system  of  equations  is  accompanied  by  the  following  ODE 
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2[b  ^  +  E  ^  -_1] 

^^?  _  -2  ^  _   2   3p 

-—-  -      a,    ■ — -_  ,  Ki) 


_2 
which  is  linear  with  respect  to  a^  for  any  Y  (ratio  of  specific  heats) 

on  the  assumption  that  p   is  given  (actually,  a  Picard  iteration  is 

taken  with  respect  to  the  nonlinear  term  in  p)  .   Here: 


3ii)l  9ij;,  d\l), 

'I'l  *?  »—  — p  ,    , 

D    '    =    ijj/ All    *    ^'l>]^i^]2    ■"    ^1^22    '  ^^^ 


E        =   ipiAip    +   v,App    ,  (6) 


li^l    is  the  toroidal  flux  at  plasma  edge  and  is   the  normalization  factor 


i|;i(v,t)    =   i|;i(t)i|;i  (v.t)    ,  ();■,  (v,t)  6:  [O.l] 

Aj^-j,  i,j  =  l,2  are  elements  of  the  susceptance  matrix,  A, 
p  is  the  pressure  which  satisfies  the  adiabatic  law 


p  =  p^e°°   =  (aiCi)^e°°  .  (7) 


p  is  the  plasma  density, 

a^  is  the  auxiliary  variable 


a,  =  _i  ,  (8) 

'    3V 

V  is  the  volume. 

Besides  independent  and  dependent  variables  and  their  derivatives  the 
function  F  (1)  depends  on  quasi  geometrical  coefficients  R^ ,  S^ ,  det  A, 
and  R, ,  see  [1,2,3]. 

2.  The  simplest  time  marching  scheme  is  the  following  explicit 
scheme  for  the  determination  of  v: 

yt+At  =  v^  +  AtH(v^)  (9) 

Here     At      is     the  time  step,    and  H(v    )    is   a  finite  difference  form  of  F 
(1).       For       a       linear,     constant    coefficient    problem      without 
lower  order  terms  the  time  step,    At,    is   determined  by 


At(t)   =    (A4^i)2  ^i^ ,  CFL   <    .5  (10) 

max   X(t) 


This  formula  is  used  in  the  current  problem.  Here: 


max  X(t)  =  max  max  X^(i,t)  ,    I   =  1,2,3  ;  i  =  l ,2 ni  (mesh  points)  ,   (11) 

i        i 

X^   are  the  eigenvalues  of  the  diffusion  matrix  A_  ,  [1,2,3],  associated 
with  second  space  derivatives  of  v  according  to 


-iln 


3v      a^v^  9^v  ^_  3v 

f(4;-,v,   -f.,  — IJ   =   A       — 1  +  G[i|iT,y,   _f-j 
3(1^^      9tl;f  Vl    3,];^  3^^ 


(12) 


Matrix  A       is  block  triangular 
^1 


0  0 


^21        ^22        ^23 
^31        ^32        ^33 


(13) 


Its  eigenvalues  are 


^    t3/2 


-2 


af  det  A  , 


(14) 


r,   ,   =  I      ^P<^'>       lR°   .   S°   ±    [(R^S°)2   -   4     "P   ^   ^^^      R?S?]l/2j      (^5) 


Yp   +    <B^> 


2    N      11 
Y(p   +    <B'^>) 


Here 


R?  = 


1        -p3/2     ^  1 


(16) 


(T-I)p(p  +  <B^>) 
Tp<B^> 


^l^k^l 
t1/2      ' 


(17) 


Er,     =      1       +     —      . 


.2   +    .8  -1    . 
^1 


(18) 
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T  is  the  plasma  temperature, 

n°  is  the  constant  corresponding  to  the  parallel 

component  of  plasma  resistivity  tensor, 
k-  is  the  absolute  heat  conductivity, 
<B  >  is  the  surface  average  of  the  magnetic  field  (squared). 

For  solving  transport  equations  (1)   the  original  code  utilizes  the 
explicit  scheme  (9)  with  the  Courant-Friedrichs-Lewy  constant 

CFL  =  .4 

This  slight  reduction  from  .5  is  made  to  provide  a  margin  of  safety 
from  the  linear  bound  (10). 

3,  Now  we  introduce  the  following  Crank-Nicolson  type  of  implicit 
scheme 


t^^+At     t^ 


For  linear,   constant  coefficient  problems,   this  scheme  should  be 

unconditionally  stable.   In  practice,  (19)  describes  a  highly  nonlinear 

t^^+At 
set  of  equations  for  v      which  are  solved  by  direct  iterations 


t^+At    t^^ 
vjl^r  =  x''  -  AtH[^J /  -   ]  (20) 


This   generally  makes  the  scheme  only  conditionally  stable.   In 
addition,  since  each  iteration  costs  as  much  as  one  step  of  the 
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expliclt  scheme,  it  is  necessary  that  the  effective  increase  of  the  CFL 
number  be  greater  than  the  number  of  iterations.  We  point  out  that  a 
similar  iteration  can  be  attempted  for  the  backward  Euler  scheme 


t,,+At    t„      ,  t,,+At 


^  -  V  "^  +  AtH[v  ^   "")  (21) 


which  is  only  first  order  accurate  in  time  but  more  dissipative.   The 
criterion  for  convergence  is 

My||2  <  e  .  (22) 


where 


/ 


ni 


1        [Vj.T^i  -  Vj.    j 


ylls  =       ^^\  ^.  I   .  (23) 


i=1 


and  G  is  specified.  However  this  condition  is  not  sufficient  and  has 
to  be  accompanied  by  some  additional  conditions.  In  order  to  describe 
them  let  us  consider  a  case  when  plasma  properties  at  the  edge  satisfy 
the  Neumann  condition,  implemented  numerically  by  zeroth  order 
extrapolation 

T(ni,t)  =  T(ni-1 ,t)  (24) 

In  this  case  during  plasma  transport  the  temperature  at  the  edge  may 
rapidly  increase.   For  the  initial  condition  of  T  =  2'10~^   and 
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E,^    =  2.10    at  the  plasma  edge,  time  interval  of  .1  and  ni  =  25, 
T(ni,t)  is  shown  in  Fig.   1.   For  this  problem 

max  A(t)  =  A^ (ni , t) 

or,  in  accordance  with  (l^l), 

max  let)  =  [n°  a^Cni.t)  det  A(ni)  JT'^-^^^^^  _  ^)        (25) 

This  curve  is  given  in  Fig.  2.  We  conclude  that  during  the  first 
.01  -  .02  units  of  time  max  A(t)  decreases  very  rapidly.  Therefore, 
the  explicit  time  step.  At  (10),  increases  rapidly.  A  smaller  time 
step  is  necessary  to  resolve  the  transient.  We  have  tried  to 
adaptively  adjust  the  time  step  according  to  the  behavior  of  the 
iteration  by  introducing  the  following  two  conditions: 

(i)     if         e^  >  My||2  >  e  .  (26) 


then 


CFL^MP  =  a  CFL™P  ,  (27) 


TMP 

where  a  <  1   is  a  specified  factor,  and  CFLq   is  an  initial  value  of 
CFlIMP; 


(ii)      if  Myll2  >  ^1  >  ^  .  (28) 


it  is  necessary  to  return  to  a  previous  time  t  =  tj^  -  At,   because  At 
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IMP 
can  exceed  its  critical  value.  Proper  values  of  e,  e-]  ,  a  and  CFL^ 

together  with  the  following  precaution  permit  us  to  avoid  instability. 

During  the  initial  transient  behavior^  where  the  solution  changes  very 

rapidly  in  time,  it  can  be  more  efficient  to  keep  a  fixed  small  CFL, 

For  the  results  below,  the  first  few  time  steps  were  run  with  CFL  =  .4. 

H.  The   proposed   implicit   scheme   includes  the  following 

parameters: 

CFLq"^^  =  the  initial  CFL, 

a  <  1  =  a  relaxation  factor, 

e,e,  >  e  =  permissible  errors  in  ||v||2. 

IMP 

n^  =  maximum  number  of  iterations  before  CFL    has  to  be  changed 

The  optimal  values  of  these  parameters  (with  regard  to  minimum  total 
iterations)  in  a  case  of  the  Neumann  condition  at  the  plasma  edge  are 


CFL™P  =  1  ,     a  =  .97  ,     "^  =  7  ,     e  =  10  ^  ,     e^  =  10  ^ 


For  this  set  of  parameters  and  time  interval  of  .1  unit  of  time,  the  real 

characteristics  of  the  implicit  scheme  are  given  in  the  following  table 
for  ni  =  25  (uniform  mesh  points)  with  initial  T(ni)  =  2*10  , 
l^  (ni)  =  2.10"^. 
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TABLE  1 


# 

TIME 

NUMBER 

OF  TIME 

STEPS 

NUMBER  OF 
ITERATIONS 

ATONE 
TIME  STEP 

CplIMP 

V 

2 

1 

0   10-5 

1  -500 

2 

.4 

109 

2 

10'^  -  .001 

501  -4-102 

2 

.7837 

10"^ 

3 

.001 -.011 

4103    8-10^ 

2 

.6730 

10"9 

4 

.011  -  .100 

8-103     1510^ 
(14667) 

2 

.5778 

10"9 

CFL^^P  as  a  time  function  and  as  a  function  in  a  number  of  time  steps 
is  shown  in  Fig.  3a  and  b,  respectively.  Since  an  average  number  of 
iterations  is  two  (Table  1).  the  implicit  scheme  can  be  preferable 
(with  regard  to  the  explicit)   if  At™P(t)  >  2AtE^P(t).  Functions 


At™P(t)  and  At^^^(t)  are  given  in  Fig, 


In  fact, 


At™P  ~  i.eAt^^P 


other  characteristics  of  both  schemes  are  presented  in  Table  2. 
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TABLE  2 


TOTAL 

TOTAL 

CPU 

NUMBER 

NUMBER 

TIME 

SCHEME 

OF  TIME 

OF 

IN 

STEPS 

ITERATIONS 

SEC 

IMPLICIT 

14647 

29576 

324 

EXPLICIT 

24144 

24144 

211 

Clearly,  the  implicit  scheme  is  ineffective. 

5.  Consider  the  case  where  the  proposed  scheme  is  much  more 
effective  than  the  explicit:  plasma  transport  with  the  Dirichlet 
condition  on  the  edge.  Besides  other  characteristics  the  temperature 
at  the  plasma  edge  is  fixed  in  this  case.  If  the  specified  edge 
temperature  is  small  enough,  the  maximum  eigenvalue  occurs  at  the  edge, 
and  the  maximum  allowed  time  step  is  determined  by  the  behavior  near 
the  plasma  edge.  In  practice,  the  explicit  scheme  will  run  somewhat 
beyond  its  stability  bound  because  the  maximum  eigenvalue  decreases 
rapidly  away  from  the  edge.  Since  temperature  at  the  plasma  edge  is 
fixed  and  therefore  time  independent, 


max  X   =   A, (ni)  , 


(29) 


TMP        TMP 

and  CFL    and  At    also  become  time  independent.   Preserving  the  same 
optimal  values  for  e,  e-,  ,  and  n^  as  above,  it  is  natural  to  increase  a: 


■1 1- 


.99  ,    e  =  10 


-ij 


c^  =  10 


"o  =  7 


In  accordance  with  (10),  (25)  and  (29)  the  optimal  CFL^"''  depends  on 
the  temperature  at  the  plasma  edge  and  on  the  number  of  mesh  points* 
and  their  regularity.  Let  the  mesh  points  be  equidistant  and  ni  =  25, 
as  above.  Characteristics  of  the  proposed  implicit  scheme  as  functions 
of  the  temperature  at  the  plasma  edge  are  given  in  Table  3. 


TABLE  3 


= 

TEMPERATURE  T 

AND  RELATIVE 

PLASMA 

DENSITY  £  , 

AT  THE 

PLASMA  EDGE 

MAX  X 

IMPLICIT  SCHEME 

EXPLICIT  SHCEME 

efficiency 
cfl'mp 

»7CFL= r— 

2CFlEXP 

CFL 

At 

NUMBER  OF 

ITERATIONS 

AT  EACH 

TIME  STEP 

^     2 

CFL 

At 

1 

1    10  ■• 

1.10    10^ 

6.1 

9.66    10"' 

2 

10^    10^ 

.4 

6.33    10^ 

7.6 

2 

5     10  2 

3.10    10"* 

10.1 

5.60    10 '' 

2 

10  ^    10^ 

.4 

2.22     10^ 

126 

3 

2    10  2 

1.22    10^ 

23.0 

3.26    10  7 

2 

10  9    10 8 

.4 

5.67     10^ 

287 

4 

1     10  2 

3.46    10^ 

36.0 

1.80    10  ' 

2 

10^    10  8 

.4 

2.00    10^ 

45.5 

5 

5    10-3 

9.80    10^ 

56.0 

9.92    10-8 

2 

10^    10  8 

.4 

7.08    10^° 

70.9 

6 

2  .  10-3 

3.87     10^ 

71.6 

3.07    10^ 

2 

10  9     10^ 

.4 

2.07     10  ''° 

89.5 

7 

1     10  3 

1.10    10^ 

78.0 

1.24    10^ 

2 

10'^    10  8 

.4 

6.36    10^  ■" 

977 

*With  ni  =  ^9   and  T  =  2-10  ^  at  the  plasma  edge  the  optimal 


IMP 


IMP 


GFL    =11.0  while  with  ni  =  25  and  the  same  T,   CFL"^"^  =  23.0. 
However,  a  systematic  investigation  of  this  was  not  made. 
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IMP 


Clearly  Hqpl   '   ■'•25CFL  is   a   prior   upper  bound  to     a     real     posterior 


'CFL 


pyp      IMP 

QpyCAr/Qpyxwr  Qj_^QQ     gach  iteration  of  the  implicit  scheme 


requires  a  little  bit  more  computations  than  one  time  step  of  the 
explicit  scheme.  As  seen  from  Table  3,  the  standard  explicit  scheme  in 

fact  cannot  be  used  when  temperature  at  the  plasma  edge  £  2«10 

—9 
because  At  becomes  of  order  10   .   Fortunately  it  was  found  that  with  the 

Dirichlet  boundary  condition  the  explicit  scheme  is  stable  with  large 
increase  in  CFL  --  Table  4. 


TABLE  4 


# 

TEMPERATURE T 

AND  RELATIVE 

PLASMA 

DENSITY  ?, 

AT  THE 

PLASMA  EDGE 

IMPLICIT  SCHEME 

EXPLICIT  SCHEME 

EFFICIENCY 

CFL  "VP 
Tlrn  - 

CFL 

At 

CFL 

At 

2CFL''^'^ 

1 

1  •  10- '' 

6.1 

9.66  -  10-7 

2.0 

3.17-10-7 

1.5 

2 

5-  10"2 

10.1 

5.60  .  10-7 

2.4 

1.34  •  10'^ 

2.1 

3 

2  •  10'2 

23.0 

3.26  •  10-7 

3.6 

5.10-  10"^ 

3.2 

4 

1  •  10'2 

36.0 

1.80-  10'^ 

3.8 

1.90  •  10"^ 

4.7 

5 

5  •  10'^ 

56.0 

9.92  ■  10'^ 

3.8 

6.73  •  10"^ 

7.4 

6 

2    10-3 

71.6 

3.07  •  10'^ 

3.8 

1.63-  10'^ 

9.4 

7 

1  •  10-3 

78.0 

1.24-  10"^ 

3.8 

6.02- 10- ""O 

10.3 

However,  as  it  follows  from  this  Table,  the  implicit  scheme  is  still  much 

TMP  FY  P 

more  efficient.   CFL     and  CFL    of  Table  ^   are  also  shown  in  Fig. 
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6.  The  accuracy  of  the  proposed  implicit  scheme  was  analyzed  by  a 
comparison  of  temperature,  pressure,  and  density  profiles  obtained 

after  .01  unit  of  time  of  plasma  transport  by 

EXP 

a)  the  explicit  scheme  with  CFL    =  3.6  and  ni  =  25  (reference 

case) , 

TM  P 

b)  the  implicit  scheme  with  CFL    =  23  and  ni  =  25,  and 

IMP 

c)  the  implicit  scheme  with  CFL    =  11  and  ni  =  49. 

Initial  data:  T^^  =  2«10~^  and  g^  ^^  =  2*10"^  fixed  in  time  (the  Dirichlet 
conditions).  Profiles  T(i(;^),  p(4j^),  and  p(ilJi)  are  presented  in  Table 
5.  As  seen  from  this  table,  discrepancies  in  corresponding  profiles  do 
not  exceed  ^%  on  the  interval  i&[0,ni-6],  i.e.,  on  "ij;^  ^  [0,  .75]  (for 
all  three  runs).  On  [ni-2,ni-l]  all  those  profiles  lose  their 
smoothness  and  have  an  artificial  "zig-zag"  (see  also  Figs.  6  and  7). 
We  guess  that  the  nature  of  those  "zig-zags"  are  due  to  a  boundary 
layer  caused  by  nearly  singular  boundary  condition.  This  models  the 
behavior  with  a  vacuum  and  a  strong  sink.  Therefore  at  the  plasma  edge 
zone  there  may  not  be  sufficient  resolution  of  the  solution. 
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Profiles   T(i|)j^),    ?{^^) ,    qO)^)        ni   =   25   and   49 


TABLE  B 


p 

1 

TEMPERATURE 

PRESSURE 

DENSITY 

EXP  25 

IMP  26 

IMP  49 

^ih 

J 

i- 

EXP  25 

IMP  25 

IMP  49 

^ih 

^{^ 

EXP  25 

IMP  25 

IMP  49 

J 

ih 

^(ih 

1 

0 

4  409 

4  406 

4  403 

.07 

14 

4745 

4741 

4736 

.08 

.19 

1076 

1076 

1076 

0 

0 

2 

02083 

4  317 

4555 

1055 

3 

04167 

4232 

4  228 

4  232 

.09 

0 

4377 

4373 

4377 

.09 

0 

1034 

1034 

1034 

0 

0 

4 

0625 

4  145 

4199 

1013 

5 

08333 

4056 

4  051 

4056 

.■"> 

02 

4022 

4018 

4023 

.10 

-.02 

09918 

09917 

09918 

0 

0 

6 

1042 

3  967 

3848 

09701 

7 

1250 

3  876 

3871 

3  876 

.13 

0 

3675 

3670 

3675 

.14 

0 

09482 

09481 

09482 

0 

0 

8 

1458 

3  786 

3505 

09261 

9 

1667 

3  694 

3689 

3  693 

.14 

03 

3339 

3334 

3338 

.15 

.03 

09039 

09038 

09039 

0 

0 

10 

1875 

3  602 

3175 

08817 

1 

2083 

3511 

3  506 

3610 

.14 

03 

3018 

3013 

3017 

.17 

.03 

08596 

08594 

08595 

02 

0 

2 

2292 

3419 

2863 

08374 

3 

2500 

3330 

3  324 

3  328 

,14 

06 

2715 

2710 

2713 

.18 

.07 

08154 

08153 

08153 

0 

0 

4 

2708 

3  237 

2568 

07933 

5 

2917 

3  149 

3  143 

3  147 

.19 

0 

2429 

2424 
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Thus,  we  conclude  that  the  implicit  scheme  possesses  the  same  level  of 
accuracy  as  the  explicit  scheme  (with  the  same  mesh  points). 

7.  Temperature  profiles  versus  geometry  steps  are  given  in  Fig. 
6  (ni  =25).  They  prove  the  statement  that  three  geometry  steps  are 
sufficient  for  convergence.  Temperature  profiles  versus  time  are  shown 
in  Fig.   7  for  t  =  0,  .025,  .050  units  of  time. 

8.  Conclusion: 

a)  the  proposed  implicit  scheme  possesses  the  same  accuracy  as  the 
explicit  scheme; 

b)  this  scheme   is  much  more  effective  in  the  case  of  plasma 

transport  with  substantial  heat  loss  (the  Dirichlet  conditions  at  the 

Fy  p 
plasma  edge)   even  if  CFL    >>  .4.  Its  efficiency  rapidly  increases 

(with  regard  to  the  explicit  scheme)  with  the  increase  of  a  singularity 

at  the  plasma  edge; 

c)  the  implicit  scheme  is  ineffective  in  the  case  of  plasma 
transport  with  the  Neumann  conditions  at  the  edge; 

d)  both  numerical  schemes  have  to  be  viewed  as  complementary  to 
each  other.  The  effective  code  must  contain  both  and  use  one  or  the 
other  depending  on  which  boundary  conditions  at  plasma  edge  are  chosen 
for  a  particular  plasma  transport  problem. 
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TIME  STEP  JvtVS  TIME 

THE  NEUMANN  CONDITION  AT  THE  EDGE 

(.1  UNITOF  TIME) 
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TEMPERATURE  PROFILES 

THE  DIRICHLET  CONDITION  AT  THE  EDGE 

(.05  UNIT  OF  TIME) 
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TEMPERATURE  PROFILES 

THE  DIRICHLET  CONDITION  AT  THE  EDGE 

(.05  UNIT  OF  TIME) 
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ON  THE  ODE  WHICH  ACCOMPANIES  A  SET  OF  PLASMA  TRANSPORT  EQUATIONS 

S.  Dinkevich* 

1 .   The  ODE  which  accompanies  a  set  of  plasma  transport  equations 


3v  3v   9^v 


9t 


follows  from  conservation  of  momentum 


:.2y    9a,  9v 

+   ai(l)o  — 

9ilj^  a^      94^^  1    ^  9ii;^ 


-  — ^  = +   a,(l)5  - —  +   afB        ,  (2) 


l'^2  ^TT-      °r 


where 


4;  9A,  ,  9A,p  9A22 

9i|j^  '      9>Jj^  ^      9(1;^ 


If  we  introduce  the  adiabatic  equation 


p  =  e°°p'^  ^   e^'oCa^Ci)^  ,  W 


we  obtain  the  following  ODE  for  computation  of  the  auxiliary- 
variable  a,  =  9i|;^/9V: 


*EBASCO  Services  Inc.,  New  York,  New  York, 
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-^— 3—  *   a^<{)2  __  +  -1:  __  +  p  __  +  afB  1  =  0     5 


With  arbitrary  Y  this  is  a  nonlinear  ODE.  Therefore  in  early  runs 
of  the  code,  Y  =  2  is  taken,  which  makes  it  linear. 

Clearly,  all  variables  which  are  contained   in  (5)  completely 

determine  pressure  p  (U).  Therefore  it  is  expedient  not  to  replace 

dp/d^^  by  \i)^  derivatives  of  a^  ,  E,^  and  o^ .  The  original  eq.  (2) 
is  already  an  ODE  for  a^  determination.  It  is  linear  with  respect 
to  a|  and  does  not  contain  Y,  i.e.,  linear  for  any  Y: 


o    2 

2(B*We*1    '^M 

9a^ 

'*-      . 

2       3p 

84-1 

' 

D^l    ^'^1 

(6) 


where 


D  ^  =  All  +  2viAi5  +  vfA^,  (7) 


Ml  ^  ^"1"12  "  ^1"22 
M2  *  ^1^22 


*1 
E   =  Ai ,  +  Vi Ao,  (8) 


3.  With  regard  to  normalized  variables 


-    *1       -     »        -     *        -    ^1 


one  can  write 
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2[b    ^    +   E    '    -— -J 
^4   _   _-2   ^  _      2       3p 


(10) 


where 


ijj^  ^p  9t P  ■if  P      Tl 

D    '    =    t|>i    A-,-,    +    2\|j^v^A^2    ■"    ^l'^22    "    '''l    ^         '  ^^^^ 

^.   In  numerical   computation  dp/d^)-^      is  most  sensitive   quantity, 
Therefore  the  computation  of  a^  requires  2-3  iterations  in 

p(k.1)  ^  /O^-JX^^JY  ^       k=0.1.2  (U) 


especially  for  first  couple  of  hundreds  of  time  steps.  Note  also 
a,  ^  >  0,  while  a,  p  ^  '^ '  negative  a^  is  meaningless  and  should  not 
be  computed  at  all. 
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ANISOTROPIC  PLASMA  TRANSPORT  IN  3-D: 

FLOW  CHARTS  OF  THE  CODE 

S.  Dinkevich* 

This  report  presents  flow  charts  and  block  diagrams  for  various 
versions  of  numerical  codes  constructed  on  the  basis  of  Grad's 
Alternating  Dimension  formulation^"^  for  3-D  anisotropic  plasma 
transport.  Detailed  formulation  has  been  reported  in  another  report.  ^ 
Notations  followed  in  this  report  are  as  introduced  in  Ref.  3. 

The  first  appearance  of  the  original  code,  which  these  versions 
are  based  on,  were  written  by  Alvin  Bayliss.  Later  a  reference 
version,  "Ref.  3-D",  was  prepared  by  John  Hogan,  Starting  from  Hogan's 
version  we  have  made  corrections,  modifications  and  additions. 

The  equilibrium  solver  that  has  been  coupled  to  the  transport 
originated  in  Ref.  5. 

Presently  we  have  two  different  versions,  "SOL  16"  and  "SOL  17". 
SOL  16  uses  implicit  time  marching  while  SOL  17  has  explicit  marching 
as  was  in  Ref.   3-D. 


*EBASCO  Services  Inc.,  New  York,  New  York. 
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FIGURE  1 
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DO    777    JJ  =  1,  JGTSTP 

JSTEP  =  JJ 
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RESPECT  TO  a] 

2)  COMPUTE  a,  =  5,,   >  0 

3)  S,    IS  RECOMPUTED  Nl  10  <  3  TIMES 

4)  RECOMPUTE  p  =  e°  (a,  J]''  AND  REPEAT  a, 
p  IS  RECOMPUTED  NI  9  <  3  TIMES 


T 


FIGURE  3-2 
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i 


callCOMPUL 


'11 


12 


22 


9A^2  9A,,        .    ij/ 


3A,^  3A,,        ,    i!;. 


8A^^  (3A  ^     ^ 

A,2^  -  ^2y^7-'^,<^ 


A^^  =  a,^(c.-f,V^;^,ff;^?C: 


4',-  A^^ 


1,  NI 


T 


CALL 

DIAG 

■^1 

=  a,  .C/;  A„+i^,   A,^) 

s 

^7 

=  S,  (vj  A,2  +  i^,  A^) 

P 

=  °J, 

P 

T 

=  e°OpT-^ 

^ 

.^0,^3. 

i  =  1,  NI 

I 


CALL  DERVYH 

30,      302      9^,       31      3p 
3C^,  '  3;/;,     d^^       3^,     3^, 

H 

i  =1,  NI 

301      302       3p 

3^,      ^^1      3i//, 

s 

1  =1,  NI 

T 


FIGURE  3-3 
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i 


CALL  JSQBSQ 

30,  30^ 


<J-B>  =  a,    0,  -^  -  0,  -^ 


<J,^>    =  <J2>  -  <J> 

It  1 


"^2    V    <J-B>^ 


<B'^>/       <B2> 


I 


CALL  EVAL 


DIFFUSE  MATRIX  A- 


V,       i  =  I,  Nl 


At    =    (A^,) 


i  =  I,  Nl 

2  CFL10 

MAX(X,   ,  X^) 


i  =  I,  Nl 


CALL  PLEDGE 

AN  AUXILIARY  SUBROUTINE  TO  CHECK 
IN  THE  TRANSPORT  EOS 

T3D1  =  T3D  -  DT  ;  TT-  =  TT  •  DRBOO^ 

i//*  =  TT-  (1  +a,T3D1  +  CLj  (T3D1)^) 


3^ 
3t 


at 


=  TT-  (1  +  a,  T3D  +  OL^  (T3D)2) 


''l.NEW 


T 


CFLI0 


.4  IF  ITDIF  <  500 
CFL  OTHERWISE 


PEDGE    =    v!', 
SIEDGE  =    i;-; 


FIGURE  3-4 
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TO  701  (FIG  3-2) 


1 


BEGIN  PLASMA  TRANSPORT  AT  FINE  POINTS  i=2,  Nl-1 


5       „ 

s 

i',  (t  +  At)   =  :.,  (t)    +     J^   c^  Tj^ 

I  (t  +  At)  =  i;  (t)  +    1^  cj  F^ 

6           r, 

a,(t  +  At)=   a,(t)+     J^   c^a,^ 

i  =  2,  NI-1 

AT  MAGNETIC  AXIS: 

LINEAR  EXTRAPOLATION,  FOR  EXAMPLE, 

OqID  =  20q(2)  -  Oq(3) 

AT  THE  PLASMA  EDGE: 

?,(NI)  =if"-°(NI)  •  /3 

«BDN'                              IFIZBDN  =  0 

?,(NI)  = 

F,(Nl-1)a,{NI-l) 

,  IF  IZBDN  =  1 

a,(Nl) 

\oN-                             IFITBDN  =  0 

T(NI)   = 

T(NI-l),                        IF  ITBDN  =  1 

2T(NI-1)  -  T(NI-2),IF  ITBDN  =  2 

T(NI) 
ao(NI)=log                   , 

(P(NI))^  ^ 

T 


^  = 

-  y.(i7,(t)  +  i7;^^) 

s 

I  - 

=    VAlM^^^'"^) 

%  '- 

y.(o,(t)  +  Oo^^'") 

i  =  1,  NI 

JCITER  >  0 


FROM  FIG  3-6 


GO  TO  A  NEW 
ITERATION 


2nd,  3rd  ... 
ITERATION 


STORE  PROFILES 

OF  PREVIOUS 
ITERATION 


IMPLICIT 
SCHEME 


IMPLICIT 
SCHEME 


FIGURE  3-5 
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TO  55  (FIG  3-5) 


JCITER  <NGRAD 


< 


1 


NORMS 


EACH  ||...|lj<e 


ANY  ||...lL>e, 


58ir 


STEADY    STATE  CHARACTERISTICS 


SS  '  ^SS  '  "ss 


v^--1/7mx)^x^^^)]'^^, 


'^ss      At  V   'o 


i 


ITDIF 

=  ITDIF  +  1 

TMPR 

=  TNPR  +  1 

T=  T 

+  DT 

JZ0  = 

MOD  (ITDIF, 

500) 

f 


IMPLICIT 
SCHEME 


NGRAD  =  INITIAL 

NUMBER  OF 
ITERATION 


MGRAD  =  MAXIMUM 

NUMBER  OF 
ITERATION 


FIGURE  3-6 
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TO  700  (FIG  3-2) 


FOR  EACH  JZ0  TIME  STEP  STORE 

ITDIF  ,  ^ss  ,  ^ss  •  °ss  ' 
TTT,  DT,  \\i>\\^.  m\^.  Ilollj, 
CFL,  #  OF  MAXX,  LOCATION 
(IN  i^, ),  MAXX 


PRINT  PROFILES 

MAXA,<J2>,  <Tj„  J2,   +7?i  J^l>        i=|,  Nl 
AT  THREE  TIME  MOMENTS 

1)  AFTER  500  TIME  STEPS 

»-  10-"  t) 

2)  AT  THE  MIDDLE  OF  TIME  INTERVAL 

3)  AT  THE  END  OF  TIME  INTERVAL 


TO  NEXT 
TIME  STEP  T<TT 


RETURN 
END 


FIGURES-? 
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REMARK:   With  regard  to  the  original  code,  written  by  A.  Bayliss 
and  modified  by  J.  Hogan,  this  code  (SOL  16,  implicit  and  SOL   17, 
explicit)  contains: 
a)  new  subroutines 

1  .   ODE  (for  any  Y) 

2.      PLEDGE    [an  auxiliary  subroutine  to  compute  9\ij*/3t 
terms   in   transpor   equations.      It   uses 
^^    =   ij'od    +   a^t   +   agt^)] 

b)     almost  completely  rewritten  subroutines 

1.      DIFFUSE: 

5  _  2 

(a)       v^(t+At)    =       I  C^v^;     5^(t+At)    =       I      c^^^; 
k=l  k=l 


Oo(t.At)   =      I      c'^o^ 


6 

I 
k=l 


This  permits  us  to  estimate  a  real  contribution  of  each 
term  to  the  total, 
(b)  implicit  scheme 

EVAL 

(a)  a  diffusion  matrix  A_  is  computed. 

h 

(b)  X^,A2,Xo  are  its  eigenvalues. 

This  permits  us  to  analyze  the  influence  of  linear  and 
nonlinear  terms  in  transport  equations  and  to  construct 
a  linear  approximation  of  transport  equations. 

INDUC 

(a)  four  equations:   a^  (ij;^  A^  ^ +v^  A^  p)  =  "fi 
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<t>2 


^1^2 


^1^22 


M2'    "22 


Aoo   are 


are  written  explicity. 

Four   variants    in   computation  of   A, ^ 

given: 

equations  1-3-^^  (original  variant) 

1-2-3 

1-2-4 

1-2-3-4  (least  squares  solution) 

(b)  least  squares  solution  is  utilized  in  the  code 
Flag's  change  any  other  variant  can  be  used. 

(c)  partially  changed  subroutines 

1.  COMPAR         6.   JSQBSQ 

2.  GEOMINT        7.   SETUP 

3.  COMPUL         8.   BLKPLT 

4.  COMOUT         9.   PLTSTR 

5.  DERVYH 

(d)  subroutines  that  are  taken  out 

1 .  GEOMNEXT 

2.  EDGE 

3 .  BCON 


By 
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